Extended Gibbs ensembles with flow 
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A statistical treatment of finite unbound systems in the presence of collective motions is presented 
and applied to a classical Lennard- Jones Hamiltonian, numerically simulated through molecular 
dynamics. In the ideal gas limit, the flow dynamics can be exactly re-casted into effective time- 
dependent Lagrange parameters acting on a standard Gibbs ensemble with an extra total energy 
conservation constraint. Using this same ansatz for the low density freeze-out configurations of an 
interacting expanding system, we show that the presence of flow can have a sizeable effect on the 
microstate distribution. 
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INTRODUCTION 



The thermodynamics of all isolated finite unbound sys- 
tems is characterized by an irreducible time dependence. 
In condensed matter physics, clusters dissociation in- 
duced by photo-ionization P, 0, H, 01 or charge transfer 
collisions [5,, ^] , cannot be studied without properly ac- 
counting for the time window of the experiment [l] . Go- 
ing down to the femto-scale, the thermodynamic prop- 
erties of nuclear systems can only be accessed through 
collisions In the cluster world the dynamical evap- 
oration may still be associated to the thermodynamics 
of the liquid- vapor phase transition [l], 0, Q niaking use 
of a time-dependent temperature within the concept of 
an evaporative ensemble 0, however in the nuclear 
collisions case the time scales can be so short that the 
reaction and decay channels cannot be decoupled, col- 
lective flows appear, and the statistical equipartition hy- 
pothesis breaks down . If in the Fermi energy regime 
and in the associated multi-fragmentation phase transi- 
tion these collective flows may be only a perturbation 
in the global energetics, this is not true at SIS energies 
where they are likely to influence light cluster formation 
by coalescence [Tlj. In the ultra-relativistic regime the 
ordered and disordered motions become comparable in 
magnitude |23], and collective flows are believed to play 
an essential role in the characteristics of the transition 
to the quark-gluon plasma observed in the RHIC data 
[H, m, III particular, correlations and recombina- 
tion of thermalized quarks from a collectively flowing de- 
confined quark plasma, is supposed to be the dominant 
mechanism for soft-hadrons production (l^ . [isj . 

In all these very different physical situations, the huge 
number of available channels and the general complexity 
of the systems under study clearly calls for a statistical 
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treatment. However, the irreducible time dependence of 
the process makes the definition of statistical concepts 
like statistical ensemble, temperature, pressure, etc. un- 
clear [T^ . If it is intuitively recognized that the pres- 
ence of incomplete equilibration and collective flows may 
be treated in a statistical framework introducing extra 
constraints p^ . the procedure is not necessarily unique. 

The inclusion of collective motion in the form of a 
radial or elliptic flow in equilibrium models has been 
treated by different authors [H M, lU, M, 
The most spread approach is to suppose a full decou- 
pling between intrinsic and collective motion and assume 
for the expanding system a standard Gibbs equilibrium 
in the local rest frame [l^, [2^ . The quality of this as- 
sumption obviously depends on the degrees of freedom 
and energy regime under study. Concerning heavy ion 
collisions, this assumption may be justified in the Fermi 
energy regime because of the limited energy percentage 
associated to directed motion [(I, and in the ultrarela- 
tivistic regime by the empirical success of hydrodynami- 
cal models [13] . Some attempts have however been done 
to explicitly include flow in the statistical treatment. 
Limiting ourselves to classical systems of interacting con- 
stituents treated as elementary de gree s of freedom, the 
empirical treatment of flow in ref. [l9| has been shown 
not to modify the correlation properties of the system. 
However other empirical approaches 0, [2l| predict that 
the presence of flow should lead to a violation of sta- 
tistical equilibrium weights, with a trend towards more 
unbound configurations. Experimental data in the nucle- 
onic regime suggest that different mechanisms may act in 
distinct energy regimes [13, HI] . 

In this paper we address the generic statistical mechan- 
ics problem of the definition of a statistical ensemble in 
the presence of a collective flow. We will use the example 
of a classical Lennard-Jones system [l^ to evaluate some 
chosen observables for a statistical isolated system sub- 
ject to a radial flow. Molecular dynamics simulations on 
the same system have already shown that flow enhances 
partial energy fluctuations |3(| and at the same time can 
act as a heat sink [sil, [s^l , cooling the system and thus 
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preventing it to reach high temperatures. We wiU show 
that in the statistical hmit it can also act as a heat bath, 
since the relaxation of the microcanonical constraint al- 
lows the isolated system to explore a larger configuration 
space. 



II. TIME DEPENDENT GIBBS ENSEMBLES 

In a recent paper [23 | we have shown that flow natu- 
rally appears in the statistical picture [13, |4l| as soon as 
wc introduce constraints which are not constants of mo- 
tion. Consider an isolated physical system characterized 
by a finite spatial extension {B?) at a given time to- Intro- 
ducing the density matrix D = Y.{n) l**^"^) P*-"^ (*("^|, 
the minimum biased microstate probability distribution 
p*^"^ is defined by 



exp (-Aoi?^) 5{e-H^, (1) 



where H is the Hamiltonian, Aq is a Lagrange multiplier 
constraining the finite size, and 



VFao {E) - (-^o-R') 5 (E - Hn) 



(2) 



(«) 



is the associated density of states or partition sum. The 
dynamical evolution of eq.([T]) at times i > tg is obtained 
from the Liouville equation dtD = —i/h[H,D] [2^, or 
equivalently from the time evolution of the constraint. 
In the Heisenberg representation 



-iAtH A2 



iAtH 



[Atr 



p=i 



(3) 



where At = (t - to) and the B^p) operators are defined 
by the recursive relation 



Bip) 



(4) 



The time dependence of the process can therefore be re- 
casted in terms of an (a priori infinite) number of extra 

constraints In the simplified case of a system of 

non-interacting identical particles 
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the series reduces to the two operators 



H,B 
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i—l 



1=1 



(5) 

■P^ } (6) 
(7) 



Then the exact density matrix is given at any time t > 

by 



SjE - H) 



1 

i 

2to 



with 



2Ao, 
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Pi-n + fi- Pi 



2An 



(8) 



P{t)^^{Aty , ^.(t) = ^At. (9) 
m m 

The diabatic evolution of an isolated initially constrained 
freely expanding system can then be described as a gen- 
eralized Gibbs equilibrium in the local rest frame 



6{E - H) 



exp 
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i=l 



-13 {t) 



{pi - mh (t) Ti 



2m 



(10) 



with a Hubblian factor linearly decreasing in time, h = 
At-i. 

These equations show that radial flow is a necessary 
ingredient of any statistical description of unconfined fi- 
nite systems in the presence of a continuum; on the other 
hand, if a radial fiow is observed in the experimental 
data, this formalism allows to associate the flow obser- 
vation to a distribution at a former time when flow was 
absent. This initial distribution corresponds to a static 
Gibbs equilibrium in a confining harmonic potential. In 
this case the infinite information which is a priori needed 
to follow the time evolution of the density matrix ac- 
cording to eq.([3]), reduces to the three observables r^, p^, 
f ■ p + p ■ r. Indeed these operators form a closed Lie 
algebra, and the exact evolution of D\g preserves it alge- 
braic structure. This treatment can be easily extended to 
non-isotropic flows \23\ introducing an initially deformed 
spatial distribution. 

It is easy to see that eq.® is still exact for an interact- 
ing system in the Boltzmann limit of purely local interac- 
tions. If the interactions are non-local at the initial time 
to, this simple solution is not exact any more and higher 
order operators play a role. Considering a flnite range 



we can see 



two body interaction V ~ 'Yliu' ^ ^ 

that the first order correction in time to the static prob- 
lem i?^^) is identical to the free problem eq.®, while 
already at the second order B^^) contains an additional 
term 

S(2)=V2^-Vln,,W(r,.') 



where 



. In the case of a harmonic inter- 



action the ij'p) operators only contain quadratic terms 



r 



and 



■ pii> , with pij/ = - Pi. 



In this case the time evolution can be taken into account 
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by a suitable time dependent temperature and the intro- 
duction of a radial flow. 

For any other interaction B^'^^ modifies not only the 
temperature but also the two-body interaction. As a first 
order approximation we can however still consider the 
statistical ansatz at the freeze-out time: 



hmrj. 



(3Vn - XRi S {Hn - E) 



(11) 



where /3, A, h are Lagrange parameters imposing a given 
value for the average thermal energy, mean square ra- 
dius and local collective radial momentum at freeze out 
through the associated equations of state 



(Eth) 
(Prir)) 



d(3 
d\ 
/3r dh 



(12) 

(13) 
(14) 



In heavy ion collisions, the values taken by these state 
variables are consequences of the dynamics. They can- 
not be accessed by a statistical treatment but have to 
be extracted from simulations and/or directly inferred 
from the data itself. In the following we take ea. (|lip as 
an ansatz for the statistical description of an expanding 
system and explore its properties within a classical sys - 
tem of A'' = 147 Lennard Jones particles of mass m [29| . 
We expect this ansatz to be reasonable in the case of 
loose interaction or moderate fiows appearing at times 
close to the freeze-out time, and in the case of a fast re- 
organization of the potential energy surface, leading to a 
decoupling of the relaxation time of the interaction and 
kinetic energy. The adequacy of eq. pTj) to describe the 
time dependent expansion of the system will be explored 
in a forthcoming paper [43 |. 

The statistical ensemble described by cq. pTj) is simi- 
lar to a Gibbs equilibrium in the local expanding frame, 
with two important differences with respect to the stan- 
dard scenario [H, of a complete decoupling between 
collective and thermal motion. First, the energy conser- 
vation constraint acts on the total energy, including flow. 
This allows energy exchanges between the thermal and 
the collective motion, and therefore can modify consid- 
erably the partitions weight, as we show below. Second, 
eq. (jlip contains a term cx which plays the role of an 
external pressure [2^ . This term is the combination of a 
positive (out-going) pressure due to the expansion, and 
a negative pressure term imposing a flnite system size at 
the freeze-out time. In turn this implies that the cor- 
rect ensemble for treating an open flowing system is not 
the usual (TV, V, T) or (iV, F, E) ensemble [23, 26, 34] but 
rather an "isobar" ensemble, where the system square 



radius is constrained only in average through a Lagrange 
parameter. This is an important point, since it is well 
known that different statistical ensembles are not equiv- 
alent in finite systems [H, . In particular only in such 
isobar ensemble the heat capacity is expected to be neg- 
ative [s^l at the liquid-gas phase transition 35], which is 
at the origin of an intense research in the nuclear multi- 
fragmentation field fs?]. It is generally assumed by sta- 
tistical models that fragment or hadron partitions are set 
within a characteristic volume (freeze-out volume) which 
may depend on the thermal energy, but does not depend 
on flow [H, [2^, m, [3^. In this case the presence of 
flow does not affect the canonical configuration space of 
the isobar ensemble. Then flow can modify the parti- 
tions only because of the modified particle correlations 
in phase space jTll, [l^, [2l|, , and because the micro- 
canonical constraint acting on the total energy leads to 
a non trivial coupling between thermal and collective en- 
ergy 



III. SYSTEMS IN A HARMONIC TRAP 

It is interesting to notice that eq. (fTT|) is formally iden- 
tical to a Gibbs equilibrium with an external harmonic 
potential U — A//?^jf|. The deep connection between 
an B? constraint and radial collective motion is shown 
by the fact that it is extremely difficult from a technical 
point of view to equilibrate a Lennard Jones system in 
a harmonic trap; this situation is referred to in the lit- 
erature as "the harmonic oscillator pathology" [H, • 



A. Dynamics of Lennard-Jones systems 

Figure 1 shows a single very long molecular dynamics 
run for the Lennard-Jones particles trapped in a har- 
monic oscillator. Even if the amplitude of the initial 
oscillations is damped by the inter-particle interaction, 
it is apparent from Figure 1 that collective oscillations 
persist over extremely long times and the ergodic limit 
does not seem to be attained. This situation is virtually 
independent of the chosen collective frequency and total 
energy. As we now show, this behavior can be under- 
stood from the closed algebraic structure of the p^, 
and r ■ p operators, which is preserved in classical me- 
chanics if commutators are replaced by Poisson brackets. 

Let us consider as above an initial condition given by 
eq.([T|) within the ideal gas H — Ek + fc-R^/2 or di- 
luted Boltzmann limit. If the only constraint on the size 
is given by the harmonic potential, the density matrix 
eq.([T]) is a stationary solution of the Liouville equation. 
If conversely the system is initialized to a different av- 
erage size through an extra constraint Aq 0, the sys- 
tem will evolve with the appearance of a collective fiow 

= - I]„ ^ [pn ■ r-n + Tn ' p„) as in cq.®. Con- 
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trary to the free case, the successive constraimng opera- 
tors do not vanish for any p > 1 and can be written 
as 



N 



1=1 

N 

^(2p+l) ^ -^(_l)P(2w) 



Pi 



2 2mk^ 

2p Pi ' fi ~^ l^i ' Pi 



(15) 
(16) 



with ll> = y/k/m. This gives at any time a density matrix 
with the same functional form as eq.([8]), with an effec- 
tive temperature (3, constraining field A, and collective 
radial velocity v oscillating in time. For the purpose of 
getting analytical results it is easier to consider an initial 
condition in the canonical ensemble 



1 



cxp ( -fioH ~ XqR" 



(17) 



The series eq.Q can be analytically summed up, and 
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FIG. 1: (Color online) Upper part: total kinetic energy per 
particle as a function of time (in Lennard- Jones units) for a 
Lennard- Jones system of 147 particles trappe d in a harmonic 
oscillator of spring constant k such that lo — ^kjm = O.Oltf^^ 
for a total energy per particle e = 2.0e. Lower part: zoom of 
the oscillation dynamics over a reduced time interval. Dashed 
line: potential energy associated to the collective oscillation 



the time dependent partition sum results Z^-^- — 



with 



Tr 



exp -^(t). 



2rn 
■ r -V r ■ p 



A(t)f2 



(18) 



The time dependent Lagrange parameters are given by 



/3W 


= /3o - ^ (cos2tjAt- 1) 
k 


(19) 


\{t) 


= ^ iPok + Ao (cos 2LuAt + 1)) 


(20) 




= — ^sin2ci;At. 

muj 


(21) 



Eq. p8|l can be interpreted as a Gibbs equilibrium in the 
rest frame of a breathing system. For classical particles 
the trace over single-particle microstates is a phase-space 
integral Tr[] = J d^r J d^p and the canonical parti- 
tion sum is readily evaluated 



(2/3A 



-3/2 



(22) 



This leads to the prediction for the time dependent be- 
havior of the different observables 



2m 



r- = -- 



dp 

^l og^/3A^ 



6mA 



2/3A 



3/3 



2/3A- 
Znii) 



9i> 2/3A - v'^m 

Introducing the expressions of /?, A, we get 



_{p^ 
2m 



3 l + a;(l + cos2cjAi)/2 



2/3o 
3 1 



1 + x 

- cos 2ujAt) /2 



^flow 



;{p- r) = 



2/3o 
3 



2/3ol 



1 



■ sin 2ujAt 



(23) 
(24) 
(25) 

(26) 
(27) 
(28) 



where x = ko/k measures the strength of the initial con- 
straint. It is clear from the inspection of Figure 1 that 
over the time scale of a collective oscillation the inter- 
particle interaction can be neglected, the total energy 
conservation constraint does not seem to play an impor- 
tant role, and the canonical free particles result eq. p^ 
appears fairly accurate. The kinetic energy do oscillate 
with the double of the oscillator frequency in phase op- 
position, this collective motion breaking the ergodicity of 
the dynamics. 



B. Microcanonical Thermodynamics 

In order to study the effect of flow for the freely ex- 
panding system, we have performed numerical molecu- 
lar dynamics calculations within the statistical ensemble 
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eq. pTj) without {h = 0) and with {h > 0) the contribu- 
tion of a radial coUective flow. To study the thermody- 
namical properties of the isobar ensemble characterized 
by a size constraint Aq , we have constructed the micro- 
canonical distribution by sorting a canonical ensemble 
(45i | of the equivalent system trapped in a harmonic os- 
cillator of spring constant k = 2Ao//3. The canonical dis- 
tributions are obtained by coupling the system to a ther- 
mostat with the Andersen technique [46J . In brief, the 
coupling is made by stochastic impulsive forces that act 
occasionally on randomly selected particles. After each 
collision, the selected particle is endowed with a new ve- 
locity drawn from a Maxwell-Boltzmann distribution at 
the desired temperature T. The combination of Newto- 
nian dynamics with the stochastic collisions generate a 
Markov chain in phase space, which under some general 
conditions generates the canonical distribution [46| . 

The resulting microcanonical thermodynamics is 
shown in Figure 2 for an oscillator constant lu — O.Oltg"^. 
Close to the liquid-gas transition temperature, the canon- 
ical calculations give rise to very wide energy distribu- 
tions and the different events (n) can be sorted in total 
energy bins 



1 



LJn 



< 




by a finite size at the freeze-out time [45|. Each sin- 
gle canonical sampling can therefore be used to access 
the microcanonical thermodynamics over a wide enthalpy 
region. The microcanonical temperature is evaluated 
in each enthalpy bin as T{H) = 2 < ex > /3 
where the average is taken over events belonging to the 
same bin. The normalized kinetic energy fluctuation 
Ak = N {{e\) - (eif)2) /T^ is also represented. The 
nice agreement between estimations obtained with dif- 
ferent canonical temperatures shows the quality of the 
numerical sampling. Figure 3 shows the dependence of 



(29) 



This energy can be physically interpreted as a free en- 
thalpy for the isolated unbound system characterized 



FIG. 2: (Color online) Microcanonical temperature (upper 
part) and normalized kinetic energy fluctuation (lower part) 
as a function of the total energy inside the harmonic oscilla- 
tor obtained from an energy sorting of canonical distributions 
corresponding to thermostat temperatures (3 — 3/(2e_f(-) as in- 
dicated. All quantities are expressed in Lennard- Jones units. 
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FIG. 3: (Color online) Microcanonical temperature (upper 
part), normalized kinetic energy fluctuation (medium part) 
and mean square radius (lower part) as a function of the total 
energy inside the harmonic oscillator for two difi'erent osciUa- 
tor strengths. The horizontal lines in the medium panels give 
the fluctuation expected in the canonical ensemble. 

the results on the oscillator strength. We can recognize 
for low to values, i.e. loose constraints on the system size, 
the first order liquid-gas phase transition. The transi- 
tion is signaled by_the backbending of the microcanoni- 
cal caloric curve 34| corresponding to a negative heat ca- 
pacity, and the associated abnormal kinetic energy fluc- 
tuation overcoming the canonical limit Afc = 3/2 (ssj . 
The consistency between the two independent signals is 
again a proof of the numerical quality of the microcanon- 
ical sampling. These results are in qualitative agreement 
with the ones obtained for the Lattice Gas model in the 
same ensemble [i^ . We can also notice that in the energy 
interval corresponding to the transition, the mean square 
radius shows a kink and a slope change at higher energies. 
The spatial extension of the unbound phase grows more 
rapidly with the energy, and at the coexistence point the 
two phases have similar spatial extensions. This means 
that in this model, contrary to the Lattice case [S^ , the 
two coexisting phases at the transition temperature can 
be populated even in an ensemble which strongly con- 
strains the volume of the system. In particular the two 
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characteristic signals of a first order phase transition in a 
finite system, namely bimodality in the canonical ensem- 
ble and negative heat capacity in the microcanonical one, 
can be observed even in the isochore ensemble [sO, IIH . 

For stronger size constraints (smaller average volumes) 
the caloric curve is monotonic, the microcanonical con- 
straint reduces fluctuations well below the canonical 
limit, and the mean square radius increases linearly with 
the energy. This signals a supercritical system. From 
these calculations the critical pressure can be roughly es- 
timated as LOc ~ 0.015^^7^. 



IV. MICROSTATE DISTRIBUTIONS IN AN 
EXPANDING ENSEMBLE 

To simulate the expanding ensemble eq. ()lip . a radial 
momentum = mhrur is added to each particle and 
a microcanonical sorting is imposed on the total energy 
including flow E' = Y.^{P^ + Pri)'^/{2m) + Elj. The 
Hubble factor h employed at different energies has been 
obtained from the measured collective velocity of the 
same system freely expanding in vacuum {Ejj^^ ) accord- 

. Since the addition 



free 
flow 



ing to /i2/2m < i?^ ^ 

of flow trivially increases the total energy E, such that 
E' = (Eth) + {Efiow) > E, the comparison between the 
calculations without flow at an energy E and those of the 
ensemble including flow at an energy E' have to be made 
such that the average thermal energy of both systems are 
similar (Eth) = E. 
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FIG. 4: (Color online) Potential energy (left side) and size of 
the largest cluster (right side) at two different thermal ener- 
gies corresponding to the bound phase (lower part) and close 
to the transition region (upper part). The filled histograms 
correspond to a static equilibrium h — while the empty 
ones flow was included according to eq. (|lip . All quantities 
are expressed in Lennard- Jones units. 

The results are shown in Figure 4 for the distribution 



of the potential energy and the size of the largest frag- 
ment recognized through the MST algorithm [29]. We 
can see that for all energies the presence of flow modifies 
the distributions in a sizeable way, leading to higher fluc- 
tuations. This is easy to understand from eg. pip if we 
consider that in the expansion dynamics only the total 
energy is conserved, meaning that thermal energy fluc- 
tuations can be compensated by collective energy fluc- 
tuations. In this sense the collective motion acts as a 
heat bath, leading to distributions similar to the canon- 
ical ones. In particular if the system has a total energy 
inside the coexistence region of the flrst order phase tran- 
sition (upper part of Figure 4) the exchange with the 
flow reservoir can allow the system to explore the two 
coexisting phases. These latter differ in potential energy 
Ablj ~ 0.4e but not in average spatial extension (see 
Figure 3) and can therefore be accessed in the same en- 
semble for a given value of the average freeze-out volume. 

This result implies that signals of phase transitions 
typical of the canonical ensemble such as bimodalities 
can be pertinent also in the microcanonical framework, 
if flow is accounted for in a thermodynamical consistent 
way. Then such signals may be accessed even in experi- 
mental situations where the deposited energy is strongly 
constrained. A possible experimental confirmation of this 
prediction in nuclear multifragmentation can be found in 
ref.[47i]. 

At this point a word of caution is in order. Our ansatz 
(jlip is exact only for a system of non-interacting parti- 
cles (or in the limit of local interactions). In the presence 
of strong correlations this ansatz supposes the system re- 
laxation time be small compared to the time-scale of the 
expansion. This should be fulfilled if the average collec- 
tive velocity (vp) is much smaller than the velocity as- 
sociated to the thermal motion (vth)- Within the ansatz 
(fTTj) the local equation of state for the radial momentum 
reads 



{pr{r)) « - ^Igf ^ = hmr 



(30) 



where we have neglected the effect of the energy- 
conserving (5- function in cq. (jlip in order to have an ana- 
lytical order-of-magnitude estimate. This leads to a col- 
lective velocity vp = h{R) which should be compared to 

the canonical estimate vth = \J 3/ (/3m). In the case of the 
upper part of Figure 4 we have vp/vth ~ 0.79 meaning 
that the quality of our approximation may be doubtful. 
It is however interesting to note that the bimodal shape 
of the distribution in the presence of flow persists also 
for smaller collective motions, as long as the energy fluc- 
tuations are of the order of the energy distance Ae^j 
between the two phases. 



V. CONCLUSIONS 

To conclude, in this paper we have presented an infor- 
mation theory based formalism allowing to include coUec- 
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tive motions in the statistical description of finite bound 
or unbound systems. Molecular dynamics simulations 
performed on a Lennard- Jones system suggest that even 
in the simplified approximation of non interacting par- 
ticles, the presence of flow can influence the microstate 
distribution in a sizeable way. Indeed the presence of a 
(non-conserved in time) collective energy component can 
play the role of a heat bath, allowing for extra configura- 
tional energy fluctuations in the total energy conserving 
dynamics. In particular, close to a first order phase tran- 
sition, this mechanism is seen to give rise to a characteris- 
tic bimodal behavior, similar to some recent experimental 



observations in nuclear multifragmentation. 
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